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Abstract 

Bulk  heat  llux  parameterization  is  an  increasingly  popular  technique  for  forcing  non-coupled  ocean  models.  If  sea  surface 
temperature  (SST)  from  the  model  is  colder  (wanner)  than  observed,  then  the  net  heat  flux  will  be  higher  (lower)  than  observ  ed;  thus, 
bulk  parameterizations  tend  to  keep  model  SST  close  to  observational  SST  on  long  time  scales.  However,  bulk  parameterizations 
imply  neither  strong  damping  of  SST  variability  nor  strong  relaxation  to  near-surface  (e  g.,  at  10  m)  air  temperature  (7;,)  This  is 
demonstrated  using  SST  simulations  from  a  0.72°  *  0.72°  cos(lat)  (longitude  x  latitude)  resolution  global  HYbrid  Coordinate  Ocean 
Model  (HYCOM)  that  does  not  include  assimilation  of  any  SST  data  or  explicit  relaxation  to  any  SST  climatology,  but  does  use  bulk 
heat  fluxes.  Results  are  discussed  when  climatological  wind  and  thermal  atmospheric  forcing  for  HYCOM  arc  constructed  from  three 
dilTerent  archived  numerical  weather  prediction  (NWT)  products:  (1)  the  European  Centre  for  Medtum-Range  Weather  Forecasts 
(HCMWF)  15-year  Re-Analysis  during  1979-1993  (ERA-15),  (2)  ECMWF  40-year  Re-Analysis  (ERA-40)  during  1979-2002,  and 
(3)  the  Common  Ocean  Reference  Experiment  Corrected  Normal  'tear  forcing  version  1.0  (CORE-CNY)  based  on  the  National 
Centers  for  Environmental  Prediction  (NCKP)  re-analysis  which  spans  1948  2002.  To  investigate  the  implications  of  the  bulk  heal 
flux  approach  as  relaxation  to  SST  and  7Y  HYCOM  SST  simulations  are  used  to  demonstrate  that  model  SST  errors  with  respect  to  the 
National  Oceanic  Atmospheric  Administration  (NOAA)  SST  climatology  do  not  look  like  7],  SST  fields  from  NWP  products.  Such  a 
result  is  confirmed  for  all  simulations  forced  with  ERA- 15,  ERA -40  and  CORE-CNY,  separately.  Overall,  global  averages  of  mean 
HYCOM  SST  error  tire  0.2  °C  ( 1 .5  °C),  0.4  °C  (1.7  °C)  and  0.6  °C  (2.3  °C)  with  respect  to  NOAA  SST  (NWT  TA)  climatology  when 
the  model  uses  atmospheric  forcing  from  ERA-15,  ERA-40  and  CORE-CNY,  respectively. 

<  2008  Elsevier  B.V.  All  rights  reserved. 
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I.  Introduction  and  motivation 

Ocean  general  circulation  model  (OGCM)  simula¬ 
tions  arc  generally  performed  using  prescribed  atmo¬ 
spheric  forcing  fields,  namely,  momentum  flux  (e.g., 
wind  stress)  and  scalar  forcing  (e.g.,  net  shortwave  and 
longwave  radiation  at  the  sea  surface,  air  temperature 
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and  air  mixing  ratio  at  10  m  above  the  sea  surface). 
1  hese  forcing  fields  arc  typically  obtained  from 
archived  NWP  products.  Examples  of  commonly-uscd 
NWP  products  include  the  European  Centre  for 
Medium-Range  Weather  forecasts  (ECMWF)  15-year 
Re-Analysis  (ERA- 1 5)  (Gibson  et  al.,  1999),  ECMWF 
40-year  Re-Analysis  (ERA-40)  (Kallberg  et  al..  2004), 
National  Centers  for  Environmental  Prediction  (NCEP) 
Re-Analysis  (Kanamitsu  et  al.,  2002),  and  the  Fleet 
Numerical  Meteorology  and  Oceanography  Center 
(FNMOC)  Navy  Operational  Global  Atmospheric 
Prediction  System  (NOGAPS)  (Rosmond  et  al.,  2002). 

All  NWPs  mentioned  above  use  relatively  sophisti¬ 
cated  boundary  layer  sub-models  to  ealeulate  sur¬ 
face  fluxes,  but  they  still  depend  on  SST  which  is 
usually  an  analyzed  (non-prognostic)  field.  The  SST 
used  is  now  accurate  enough  that  errors  in  surface  fluxes 
at  the  NWP  grid-scale  are  dominated  by  errors  in 
atmospheric  fields,  such  as  wind  speed  and  cloudiness. 
These  errors  can  be  large,  and  even  the  best  surface 
heat  flux  fields  from  all  sources  (NWP  products  or 
observation-based  climatologies)  do  not  give  a  closed 
global  heat  budget.  In  addition,  using  climatological 
estimates  of  total  heat  flux  to  force  an  ocean  model 
usually  results  in  unrealistic  model  SST,  presenting  an 
inconsistency  with  the  imposed  surface  flux  (Bamicr 
ct  al.,  1995). 

Even  perfect  fluxes,  with  a  closed  global  heat  budget, 
cannot  be  used  as  the  only  forcing  for  a  stand-alone 
ocean  model  because  the  model’s  “climatology'1  is  not 
perfect,  so  such  fluxes  will  lead  to  SST  drift  (c.g., 
I  lughes  and  Weaver,  1996).  The  conventional  approach 
to  correcting  this  drift  is  to  add  relaxation  to  observed 
SST,  either  in  place  of,  or  in  addition  to,  prescribed  total 
heat  fluxes  (c.g.,  Barnicr  et  al.,  1995;  Maltmd  et  al., 
1998).  As  stated  by  Kill  worth  et  al.  (2000),  researchers 
at  this  time  did  not  possess  consistent  surface  fields  with 
which  to  force  their  models.  However,  the  use  of  a  bulk 
parameterization  has  become  an  alternative  approach 
because  it  requires  no  explicit  relaxation  to  observed 
oceanic  quantities  (e.g.,  Wallcraft  et  al..  2003;  Kara 
et  al.,  2003).  In  addition,  using  bulk  paranictcrizadons 
(e.g.,  Parkinson  and  Washington,  1979),  regional 
coupled  ice-ocean  models,  such  as  those  in  the  Baltic 
Sea,  have  been  integrated  over  decades  (Lehmann  and 
Flinrichscn.  2000). 

Bulk  parameterizations  use  the  difference  7'a  SST, 
where  7a  is  the  air  temperature  and  SST  the  ocean  model 
sea  surface  temperature.  The  air  temperature  and  other 
parameters  (air  density  at  the  air-sea  interface,  mixing 
ratio,  and  wind  speed  at  1  0  m  above  the  sea  surface)  arc 
typically  from  NWP  products  and  arc  used  to  calculate 


the  exchange  coefficients  for  latent  and  sensible  heat 
fluxes  (e.g.,  Kara  et  al.,  2002). 

The  most  important  criticisms  of  bulk  parameteriza¬ 
tions  are  that  ( I )  an  infinite  heat  capacity  exists  between 
ocean  and  atmosphere,  i.e.,  it  docs  not  give  a  closed 
system,  (2)  since  the  bulk  formula  can  be  linearized  as  a 
function  of  a  difference  between  equilibrium  tempera¬ 
ture  and  SST  (Hancv,  1971;  Han,  1984;  Paiva  and 
Chassignct,  2001),  it  can  smear  the  model  SST  when 
the  model  resolution  is  finer  than  the  forcing,  and 
(3)  The  difference,  i.e.,  Ta  SST,  generally  represent 
the  boundary  layer  physics  of  the  NWP  products,  but 
may  be  a  poor  approximation  to  observations.  Note 
that  the  statement  in  (1)  typically  holds  if  the  atmo¬ 
sphere  temperature  is  kept  constant  over  a  longer  period, 
for  example  if  climatological  v  alues  are  used.  Otherwise 
the  atmosphere  temperature  exhibits  a  time  tendency 
too,  which  should  reflect  the  effect  of  the  heat  flux. 

In  this  paper,  we  examine  the  impact  of  bulk  heat  flux 
parameterization  on  the  climatological  SST  produced  by 
simulations  from  a  particular  global  OGCM,  with  no 
assimilation  of,  or  explicit  relaxation  to,  any  observed 
SST  climatology  In  particular,  reasons  for  preferring  the 
bulk  parameterization  rather  the  net  heat  flux  plus  an 
explicit  relaxation  to  SST  are  discussed. 

A  few  salient  reasons  for  using  SST  from  an  OGCM 
to  discuss  the  bulk  heat  flux  versus  the  direct  relaxa¬ 
tion  approach  arc  (!)  SST  plays  an  important  role  in 
air  sea  interactions  through  the  net  heat  flux  at  the 
sea  surface  (e.g.,  Alexander  and  Scott,  1997);  (2)  it  is 
typically  used  in  the  bulk  heat  flux  formulation  to 
parameterize  the  stability  (Kara  et  al.,  2000);  (3)  it 
plays  a  major  role  in  controlling  long  term  climate 
variations,  such  as  the  North  Atlantic  Oscillation 
(Paeth  ct  al.,  2003);  (4)  it  is  one  of  the  best  observed 
variables  of  the  upper  ocean  ov  er  the  global  ocean; 
therefore,  SSTs  simulated  from  an  OGCM  can  be 
easily  validated  over  the  globe  (SchopI  and  I  oughc, 
1995).  In  addition,  SST  from  an  atmospherical  ly- 
foreed  OGCM  (no  assimilation  of  any  ocean  data, 
including  SST)  is  used  because  one  of  our  goals  is  to 
demonstrate  that,  although  the  bulk  heat  flux  formula¬ 
tion  forcing  an  OGCM  surface  temperature  field 
includes  the  knowledge  of  air  temperature  near  the 
sea  surface,  such  use  does  not  imply  an  improper  form 
of  SST  restoring  w  ithin  the  model. 

The  paper  is  organized  as  follows.  Section  2 
discusses  the  traditional  bulk  heat  flux  approach  for 
OGCMs.  Section  3  gives  details  of  the  OGCM  used  in 
this  study.  Section  4  presents  SST  results  from  an 
OGCM  in  relation  to  the  bulk  heat  flux  parameteriza¬ 
tion,  and  Section  5  gives  conclusions  of  this  paper. 
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2.  Bulk  heat  flux  approach 

The  total  heat  flux  ((J„cl)  at  the  sea  surfaee  ean  be 
expressed  as  follows: 

£?iwt  —  (?sw  (?iw  +  Qi  t  Os  ■  { l ) 

where  (9SVN  is  the  net  shortwave  radiation  at  the  sea  surfaee, 
(7lu  is  the  net  longwave  radiation  at  the  sea  surface,  Qi 
is  the  latent  heat  flux,  and  (7S  is  the  sensible  heat  flux. 

All  these  individual  components,  and  their  total,  Qnct, 
are  typically  avai table  from  the  archived  real-time  and  re¬ 
analysis  data  sets  or  from  climatologies  (ERA-1 5,  ERA-40 
and  CORE-CNY  are  used  here).  Shortwave,  Qsw*  and 
longwave,  Q\^  radiation  can  be  directly  measured  at  the 
surfaee,  but  there  are  far  too  feu  such  observations  to 
constrain  even  regional  fluxes.  Thus,  in  practice  all 
components  are  estimated.  All  components  are  affected 
by  surface  type  i.c ,  land  versus  sea  versus  sea-ice  (Garratt 
et  al  ,  1 998).  but  here  we  will  only  consider  the  open  ocean. 
Both  shortwav  e  mid  longvvav  e  radiation  depend  heavily  on 
cloudiness,  and  all  flux  components  except  shortwave 
radiation  depend  on  SST  (c.g.,  (iill  1982.  Glccklcr  and 
Wearc,  1997).  HYCOM  uses  a  penetrating  solar  radiation 
scheme  (Kara  et  al ,  2005a)  that  accounts  for  spatial  and 
temporal  water  turbidity'  (Kara  et  al..  2005b, c) 

While  SST  drift  in  OGCMs  may  be  significantly 
reduced  by  relaxation  to  observed  SST  (c.g..  Banner 
et  al.,  1995),  such  an  approach  has  the  important  short¬ 
coming  that  the  correct  result  is  prescribed  in  the  model 
simulation.  Rather  than  using  a  direct  SST  relaxation 
scheme,  an  alternative  approach  is  to  calculate  fluxes 
from  a  bulk  parameterization  based  on  ncar-surfaec 
atmospheric  fields  using  only  n ear-surface  wind  speed 
(va).  or  also  including  air  mixing  ratio  (#.,),  near-surface 
air  temperature  {Tlx)  (all  of  which  are  10  m  above  the  sea 
surface),  mixing  ratio  for  sea  water  (</*)  and  model  SST, 
as  explained  in  Kara  et  al.  (2005d). 

Bulk  parameterizations  are  based  on  statistical  fits  to 
observations.  They  arc  usually  inexpensive  to  calculate, 
and  the  best  parameterizations  arc  generally  very 
accurate.  Examples  of  typical  bulk  parameterizations 
for  latent  and  sensible  heat  fluxes  can  be  found  in 
DeC'osmo  el  al  (1996),  Kara  et  al.  (2000,  2002),  and 
1  airall  el  al.  (2003).  The  formulation  used  in  HYCOM  is 
briefly  described  as  follows: 

Qs  C\Cppayu(Tn  -  SST)  (2) 

Qi  -  C,LP ava(*,  ch)  (3) 

The  air  density  at  the  air  -  sea  interface  (pa)  in  kg  m  * 
is  determined  using  the  ideal  gas  law'.  The  exchange 


coefficients  (C/  and  C's)  used  in  HYCOM  include  the 
effects  of  boundary  layer  stability,  and  are  expressed  as 
polynomial  functions  of  Ta  SST,  va,  and  cjA  t/s,  the 
latter  through  relative  humidity  (RH),  where  RH  is  at  the 
air-sea  interface  (sec  http://vv\vw7320. nrlssc.navy.mil 
nasccnascc.html).  In  the  polynomial  functions,  the 
effects  of  water  vapor  flux  in  calculating  the  exchange 
coefficients  are  taken  into  account  through  RII  effects 
that  are  especially  important  at  low  wind  speeds  (Kara 
et  al.,  2005d).  Ihe  total  heat  flux  (Eq.  (1))  is  therefore 
expressed  as  a  polynomial  function  of  the  model  SST 
and  can  be  linearized  to  first  order  as 

(?nc  =  /(71*  SST),  (4) 

in  which  the  apparent  equilibrium  temperature  7*  and 
the  relaxation  coefficient  /.  are  time  and  space 
dependent  and  ean  be  computed  from  the  atmospheric 
fields  (Haney,  1971;  Han.  1984;  Paiva  and  Chassignet, 
2001).  We  have  included  (9SVN  in  the  Haney  parameter¬ 
ization,  as  is  usually  done,  even  though  it  docs  not 
depend  on  SST.  The  total  heat  flux  therefore  implies  a 
temperature  relaxation,  but  not  toward  the  atmospheric 
temperature  nor  toward  a  “correct"  SST  (7‘).  The 
equilibrium  temperature  T*  indeed  must  differ  from  Tc, 
as  it  would  otherwise  imply  zero  heat  flux  when  the 
model  SST  is  equal  to  Tc. 

2.1  Preference  for  a  bulk  parameterization 

Instead  of  using  a  direct  relaxation  to  an>  SST 
climatology,  a  common  alternative  has  been  to  use  NWP 
net  flux  plus  an  explicit  relaxation  to  the  correct  SST,  If 
Since  SSI  is  one  of  the  best  observed  geophysical  fields, 
why  should  a  bulk  parameterization  be  preferred  over 
the  NWP  net  flux  plus  an  explicit  relaxation  to  SST? 
One  answer  is  that  the  c- folding  time  implicit  in  accurate 
bulk  parameterizations  is  representative  of  the  air  sea 
interface.  It  is  possible  to  construct  explicit  relaxation 
terms  that  are  patterned  on  bulk  parameterizations,  or 
alternatively  to  use  measured  climatological  c-folding 
times  from  observations  or  NWP  fields  (or  the  results  of 
a  bulk  parameterization  applied  to  NWP  fields).  These 
approaches  are  all  of  the  Haney  type  (1  laney,  1971;  Han, 
1984;  da  Silva  et  al.,  1994;  Paiva  and  Chassignet.  200 1 ) 
as  follows 

Qi  SST)  =  Q(Ts)+^\Tt(Te  SST).  (5) 

where  the  effective  SST  relaxation  c-foldmg  time  scale 
implied  by  depends  on  many  factors.  However,  it 
is  certainly  positive  and  typically  large  enough  that 
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model  SST  docs  not  exhibit  long  term  drift.  In  practice 
relaxation  terms  arc  often  simpler  than  this,  i.e.,  a 
constant  e- folding  time  and  some  approximation  to  Tc. 

Eddy-resolving,  atmospherically -forced  ocean  models 
without  ocean  data  assimilation  should  not  relax  strongly 
to  synoptic  observed  SST  because  the  observed  SST  will 
include  SST  anomalies  associated  with  eddies  and 
western  boundary  currents  at  different  locations  than  in 
the  model  simulation.  This  is  because  of  flow  instabilities 
in  the  model  which  arc  not  a  deterministic  response  to  the 
forcing.  The  bulk  parameterization  does  not  use  the  NWP 
SST  (7's)  explicitly,  and  the  NWP  fields  arc  often  on  a 
much  coarser  grid  than  the  ocean  model  one  that  barely 
resolves  such  details  as  eddy  and  ocean  current  locations. 
However,  this  raises  the  question  of  to  what  degree  the 
implied  relaxation  in  the  heat  flux  can  smear  the  model 
SST  when  the  model  resolution  is  finer  than  the  forcing.  It 
also  raises  the  question  as  to  whether  the  NWP  atmo¬ 
spheric  temperature  7"a  is  strongly  correlated  to  the  NWP 
SST  Ts,  therefore  introducing  oceanic  information  in  the 
heat  flux  formulation.  In  contrast,  eddy-resolving  stand¬ 
alone  ocean  models  that  assimilate  sea  surface  height 
(from  satellite  altimeters)  have  ocean  fronts  and  eddies  in 
the  observed  location  (Smedstad  et  al.,  2003;  Chassignet 
et  al.,  2006;  Shriver  et  al..  2007)  and  therefore  can  relax 
directly  to  high  resolution  observed  SST,  although  they 
may  instead  (or  in  addition)  directly  assimilate  SST.  From 
the  preceding  discussion,  simple  relaxation  is  never 
appropriate  as  the  only  heat  flux  term,  so  it  must  be 
augmented  either  with  the  NWP  net  heat  11  ux  or  a  bulk 
parameterization  heat  flux. 

Regions  where  the  ocean  model's  mean  SSTs  are  less 
accurate  can  often  be  traced  to  deficiencies  in  the 
atmospheric  fields.  In  particular,  systematic  biases  in 
wind  speed  will  always  lead  to  poor  fluxes  from  the  bulk 
parameterization  and  biases  in  cloudiness  will  have  a  large 
effect  on  shortwave  radiation  (partially  offset  hy  changes 
to  downward  longwave  radiation).  Biased  winds  arc 
unlikely  to  provide  better  surface  heat  fluxes  from  the 
sophisticated  boundary  layer  sub-model  in  the  NWP  than 
from  the  bulk  parameterization.  Assimilation  of  wind 
speeds  from  satellites  should  improve  the  accuracy  of 
NWP  surface  \v  inds,  but  there  are  still  large  biases  in  some 
coastal  regions,  and  there  may  be  smaller  systematic 
biases  on  larger  scales.  Cloud  cover  is  difficult  to  obtain 
accurately,  and  known  to  be  deficient  regionally  in  both 
archived  real-time  and  re-analysis  products. 

7.2.  Relationship  betw  een  air  temperature  and  SST 

In  order  to  answer  the  question  as  to  whether  the 
atmospheric  temperature,  is  strongly  correlated  to 


the  SST,  therefore  introducing  oceanic  information  into 
the  heat  flux  formulation,  the  long  term  (climatological) 
mean  of  the  difference  betw  een  near-surface  Ta  and  SST 
from  an  observation-based  data  set  and  the  SS 1  used  by 
NWP  products  are  investigated  over  the  global  ocean 
(Fig.  1).  Specifically,  the  observation-based  data  set  is 
the  Comprehensive  Ocean  Atmosphere  Data  Set 


Climatological  mean  bias:  To  7\ 


Fig.  I .  Climatological  mean  of  the  difference  between  near-surface  air 
temperature  ( 7^)  and  sea  surface  temperature  (i.c.,  Ta  /,)  over  lhc 
global  ocean  Both  Tlt  and  T%  fields  from  COADS,  FRA- 15.  FRA -40 
and  CORF-CNY  arc  interpolated  to  the  0.72°  IIYCOM  grid,  and  the 
difference  between  lhc  two  is  formed  al  each  grid  point. 
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(COADS)  climatology  based  largely  on  ship  observa¬ 
tions  and  buoy  measurements  during  1945  1989  (da 
Silva  et  al  ,  1 994)  and  the  NWP  data  products  are  the  re¬ 
analysis  products,  ERA-15  (1979-1993),  ERA-40 
( 1979-2002)  and  CORE-CNY  (1948  2002).  GOADS 
is  the  only  observation-based  climatology  among  these 
products,  and  is  intended  to  complement  the  compar¬ 
isons.  This  is  the  new  1/2°  x  1/2°  climatology  based  on 
the  atlas  of  surface  marine  data.  Supplement  B  (http:// 
www.nodc.noaa. go\  OC5  asmdnew.html).  Note  that 
hereinafter  the  GOADS  SST,  as  well  as  the  SSTs  used 
in  the  NWP  products  will  be  denoted  as  7S. 

Although  the  time  periods  over  which  each  climatol¬ 
ogy  (GOADS,  ERA- 1 5,  ERA-40  and  NGEP)  is 
constructed  are  different,  they  are  all  >  15  years  which 
is  adequate  to  represent  a  long  term  mean  over  the 
global  ocean.  Near-surface  Ta  and  Is  data  for  each 
product  arc  available  from  the  National  Gentcr  Atmo¬ 
spheric  Research  (NGAR)  web  site  (available  online  at 
http:  dss.ucar.edu  catalogs  ).  GOADS  provides 
monthly  mean  7a  7\  fields  directly,  while  for  the 
three  NWP  products  we  form  climatological  means 
using  7a  and  the  Ts  fields  archived  at  sub-daily  intervals. 
Near-surface  atmospheric  temperature  7a  from  ERA- 15, 
ERA-40  and  NGEP  is  at  2  in,  while  that  from  GOADS  is 
at  10  m  above  the  sea  surface.  The  adjustment  from  2  to 
10  m  T,  is  very  small  as  indicated  by  the  comparison  of 
7a  values  in  I  able  I  and  hence  was  ignored. 

When  sea-ice  is  present,  surface  Ta  is  over  sea-ice  in 
the  re-analysis  products,  but  not  in  GOADS.  The  Ta  field 
from  NGEP/NGAR  is  actually  obtained  from  the 
GORE-GNY  (Corrected  Normal  Year)  data  set  which 

table  1 

Global  averages  for  near- surface  air  lemperature  (T,).  Ts  and  ditYcrencc 
between  (he  iwo  from  lhe  observational-based  ('GAOS  data  sel  and 
lhrcc  NWP- -based  products 


Data 

7« 

/; 

7a  7; 

COADS 

21.4 

22  0 

0.6 

!:RA-I5 

214 

22  1 

-0  7 

HRA-40 

21.2 

22  1 

-0  9 

COARF  CNY 

20.9 

22.0 

I  1 

Further  information  about  COADS,  FRA- 1 5.  HRA-40  and  CORF.- 
C  NY  can  be  found  online  at  several  web  sites.  For  example,  as  of  this 
writing,  the  web  addresses  are  http:  mdl. Ideo. eolumhia  edit  ' 
SOURCES  .DASILVA  SMD94  for  COADS.  http:,'  hadc.nere.ac.uk 
datiicennvf-cra  I.RA.hlml  for  HRA-15,  http  Kuk.iierc.at.uk  data 
eemwt-e40c40  background,  html  for  FIR  A  -40,  and  http:  datalgtdl 
noaa.gov/nomads/fonns  mom4'COR I  CNYI  IpO.liiml  for  OORF.- 
CN  Y.  The  web  addresses  may  be  subject  to  changes  by  the  originators. 
All  the  NWP  products  have  ditTcrcnl  boundary  layer  paramctcriza- 
tions,  ph>sics,  data  assimilation  methods  and  different  satellite  data 
used  in  the  assimilations.  Therefore,  differences  in  their  output 
variables  arc  expected. 


corrected  for  excessively  cold  values  in  the  Antarctic. 
Ignoring  high  latitudes  where  sca-icc  forms,  the  fields 
arc  broadly  similar  to  each  other  with  climatological  Ts 
warmer  than  T:i  nearly  everywhere  but  usually  by 
<  I  °G.  The  exception  is  in  areas  where  strong  currents 
transport  water  which  is  significantly  wanner  than  the 
air  in  the  mean,  c.g..  Gulf  Stream,  Kuroshio  and 
Agulhas.  Further  details  can  be  found  in  Kara  et  al. 
(2007). 

The  same  long  tenn  mean  fields  shown  in  Fig.  I  arc 
plotted  as  scatter  diagrams  of  T;i  versus  Ts  oxer  the 
global  ocean  (Fig.  2).  The  regions  where  ice  is  located 
arc  excluded  in  the  evaluation  procedure.  There  is  a 
strong  linear  relationship  between  TA  and  SST  with  R 
values  >0.99  for  all  data  sets.  This  is  clearly  evident 
from  the  slope  values  being  close  to  I  in  all  cases  ( 1 .008, 
0.991,  1.012  and  1.033  for  COADS,  ERA-15,  ERA-40 
and  GORE-GNY,  respectively).  Intercept  values  are  also 
similar  with  values  of  0.468,  0.911,  0.711  and  0.405. 
The  difference  between  Ta  and  Ts  (i.e.,  T.x  Ts)  ranges 
from  -  0.64  °G  for  GOADS  to  -  1 .09  °C  for  CORE-CNY 
(  Table  I). 

In  contrast,  near-surface  Tlx  and  Tu  7!,  are  only 
weakly  correlated  (Fig.  3)  with  linear  correlation  coef¬ 
ficients  of  0. 1 2,  0. 1 4,  0. 1 7  and  0.25  for  GOADS, 
ERA- 15,  ERA-40  and  CORE-CNY.  None  of  the  cor¬ 
relation  values  are  statistically  significant  in  comparison 
to  a  correlation  value  of  0  at  a  95%  confidence  interval 
based  on  the  student’s  /-test.  Thus,  7],  itsclt  does  not 
have  strong  influence  on  TA  Fs. 

3.  HYCOM  description 

The  HYbrid  Coordinate  Ocean  Model  (HYCOM)  is 
based  on  a  primitive-equation  formulation  discussed  by 
Bleck  (2002),  in  detail.  There  have  been  several 
HYCOM  applications  investigating  a  variety  of  pro¬ 
cesses  in  different  ocean  basins  and  enclosed  seas. 
Examples  of  such  studies  include  the  North  Atlantic 
(Chassignct  ct  al..  2003;  Flalliwcll  2004;  Thacker  ct  al., 
2004),  the  Indian  Ocean  (Han  ct  al  2004),  the  tropical 
Pacific  (Sliaji  et  al.,  2005),  and  the  Black  Sea  (Kara 
et  al.,  2005a,b,c).  HYCOM  has  also  been  used  as  the 
ocean  component  of  a  coupled  atmosphere  ocean 
model  in  global  climate  studies  (e.g..  Sun  and  Hansen 
2003:  Bleck  and  Sun  2004),  and  for  eddy-fresolving 
ocean  prediction  (Chassignct  et  al.,  2006). 

For  the  present  study,  we  will  introduce  HYCOM 
configured  for  the  global  ocean,  including  the  latest 
advances  in  the  model  development  (sec  Appendix  for 
details).  The  model  presented  here  is  a  stand-alone 
ocean  model  with  no  assimilation  of  any  ocean  data. 
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including  SST,  and  no  relaxation  to  any  other  data 
exeept  sea  surface  salinity  (SSS)  to  keep  the  evaporation 
minus  precipitation  balance  on  track.  General  features  of 
the  global  atmospherically-forced  HYCOM  are  given  in 


T  (°C) 

a 


u 

s 


4 

2 

0 

-2 

-4 


i  1  i  t  ^~r 
;  CORE-CNY 


I  i  i . i . iii 


0  10  20  30 


Fig.  2  Scatter  plots  of  TA  versus  Ts  based  on  annual  mean  values,  each 
at  1°  bins  over  the  global  ocean.  Both  7*  and  T%  fields  arc  from 
C OADS,  HRA-1 5.  KRA-40  and  t  ORL-CNY,  respectively 


1  ig.  .“V  Scatter  plots  of  ncar-surfacc  air  temperature  (7a)  versus  7a 
based  on  annual  mean  values  (see  Fig  I),  each  at  1°  bins  over 
global  ocean 
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Section  3.1  and  the  atmospheric  forcing  used  in  the 
model  simulations  is  explained  in  Section  3.2. 

3.1.  General  features  of  global  HYCOM 

HYCOM  contains  five  prognostic  equations:  two  for 
the  horizontal  velocity  components,  a  mass  continuity  or 
layer  thickness  tendency  equation  and  two  conservation 
equations  for  a  pair  of  thermodynamic  variables,  such  as 
salt  and  potential  temperature  or  salt  and  potential 
density  (Block,  2002).  The  model  behaves  like  a 
conventional  a  (terrain-following)  model  in  very 
shallow  oceanic  regions,  like  a  z-levcl  (fixed-depth) 
coordinate  model  in  the  mixed  layer  or  other  unstratified 
regions,  and  like  an  isopyenie-coordinatc  model  in 
stratified  regions.  However,  the  model  is  not  limited  to 
these  coordinate  types  (Chassignet  et  ah,  2003). 
Typically,  HYCOM  has  isopycnal  coordinates  in  the 
stratified  ocean  but  uses  the  layered  continuity  equation 
to  make  a  dynamically  smooth  transition  to  z-lc\els  in 
the  unstrati  lied  surface  mixed  layer  and  a-lcvcls  in 
shallow  water.  The  optimal  coordinate  is  chosen  every 
time  step  using  a  hybrid  coordinate  generator.  Thus,  the 
model  automatically  generates  the  lighter  isopycnal 
layers  needed  for  the  pycnocline  during  summer,  but 
during  winter  the  same  layers  define  fixed-depth  z-lcvcl 
coordinates. 

HYCOM  domain  used  in  this  paper  spans  the  global 
ocean  from  78°S  to  90°N.  It  has  a  0.72°  equatorial 
Mercator  grid  between  78°S  47°N,  with  an  Arctic  bi¬ 
polar  patch  above  47°N.  Average  zonal  (longitudinal) 
grid  resolution  for  the  0.72°  global  model  varies  from 
~80  km  at  the  equator  to  ^  60  km  at  mid-latitudes  (c.g., 
at  40°  N).  I  hc  meridional  (latitudinal  )  grid  resolution  is 
doubled  near  the  equator  to  better  resolve  the  equatorial 
wavc-quide  and  halved  in  the  Antarctic  for  computa¬ 
tional  efficiency.  Hereinafter,  the  model  resolution  will 
be  referred  to  as  0.72°  for  simplicity.  Zonal  and 
meridional  array  lengths  are  500  and  457,  respectively. 
At  this  resolution,  coastal  regions  are  not  represented  in 
great  detail  (so  sigma-levels  arc  not  used),  and  the  model 
land  sea  boundary  is  at  the  50  m  isobath  (with  a  closed 
Bering  Strait). 

All  global  simulations  use  sigmaO  (sigma  theta),  i.c. 
potential  density  referenced  to  the  surface  (0  dbar)  with 
no  thennobaric  correction.  There  are  26  hybrid  layers  in 
the  vertical  In  general,  the  model  needs  fewer  vertical 
coordinate  surfaces  than,  say,  a  conventional  r-Icvcl 
model,  because  isopycnals  are  more  efficient  in 
representing  the  stratified  ocean,  as  discussed  in  Hurlburt 
et  al.  ( 1996)  and  Kara  ct  al.  (2005a).  The  target  density 
values  for  the  isopycnals  and  the  decreasing  change  in 


density  with  depth  between  isopycnal  coordinate 
surfaces  are  based  on  the  1/4°  Generalized  Digital 
Environmental  Model  (GDEM)  climatology  (NAVO- 
CEANO,  2003).  The  density  difference  values  were 
chosen  so  that  the  layers  tend  to  become  thicker  with 
increasing  depth,  with  the  lowest  abyssal  layer  being  the 
thickest.  The  near-surface  z-lcvcl  regime  is  a  natural 
consequence  of  HYCOM's  minimum  layer  thickness.  In 
this  ease,  the  minimum  thickness  of  layer  1  is  3  m,  and 
this  minimum  increases  1 . 1 25 x  per  layer  up  to  a 
maximum  at  12  m,  and  target  densities  are  chosen, 
such  that  at  least  the  top  four  layers  are  in  z-lcvel 
coordinates. 

The  simulations  use  realistic  bottom  topography 
constructed  from  the  NRL  2  min  resolution  bathymetry. 
Extensive  quality  control  of  bottom  topography  was 
performed  in  straits  and  near  coastlines.  Monthly  mean 
temperature  and  salinity  from  the  GDEM  climatology  in 
August  arc  used  to  initialize  the  model.  There  is  a 
relaxation  to  monthly  mean  SSS  from  the  Polar  science 
center  Hydrographic  Climatology  (PHC).  The  PI  1C 
climatology  is  chosen  for  its  accuracy  in  the  Arctic 
region  (Steele  ct  al.,  2001).  The  reference  mixed  layer 
thickness  for  the  SSS  relaxation  is  30  m  (30  days  in  30  m 
e-folding  time).  Ihc  actual  c-folding  time  depends  on 
the  mixed  layer  depth  (MED)  and  is  30  *  30/MLD  days, 
i.e.  it  is  more  rapid  when  the  MLD  is  shallow  and  less 
so  when  it  is  deep.  Such  a  relaxation  is  necessary  to 
prevent  SSS  drift,  and  is  in  addition  to  the  evaporation 
precipitation  budget. 

HYCOM  has  several  mixed  layer/vertical  mixing 
options  (see  Halliwcll  (2004)  for  a  discussion  and 
evaluation).  In  this  paper,  the  non-slab  K-Profile 
Parameterization  (KPP)  of  Large  et  al.  ( 1997)  is  used. 

3. 2.  A  tmospheric  forcing 

The  model  reads  in  the  following  time-varying 
atmospheric  forcing  fields:  wind  forcing  (zonal  and 
meridional  components  of  wind  stress,  wind  speed  at 
10  m  above  the  sea  surface)  and  thermal  forcing  (air 
temperature  and  air  mixing  ratio  at  10  m  above  the  sea 
surface,  precipitation,  net  shortwave  radiation  and  net 
longwave  radiation  at  the  sea  surface).  1  he  wind/ 
thermal  forcing  was  constructed  from  three  NWP 
products:  (i)  1  125°  *  1.125°  ERA- 15,  (ii)  1.125°* 
1.125°  ERA-40  and  (iii)  1 .875° *  1 .875°  CORE-CNY. 
When  using  each  product,  the  goal  is  to  include  high 
temporal  variability  in  the  forcing,  while  maintaining  a 
climatology. 

The  original  atmospheric  forcing  data  set  from  ERA- 1 5 
spans  the  period  I  January  1979  31  December  1993. 
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The  climatological  ERA- 15  forcing  used  here  consists  of 
1 5-ycar  monthly  averages  over  1 979  1 993,  interpolated  in 
time  to  6  hourly  and  with  6  hourly  sub-monthly  anomalies 
from  operational  ECMWF  in  September  1994  to  Septem¬ 
ber  1995  added  to  the  winds  only.  Choosing  another  time 
period  for  the  6  hourly  wind  anomalies  (other  than  1 994 
95)  did  not  make  any  significant  impact  on  the  model  SST. 
ECMWF  used  a  spectral  model  to  generate  the  ERA-15 
dataset.  Hie  ERA- 15  re-analysis  project  incorporated  a 
number  of  in-situ  and  satellite-based  data  types  over  the  1 5- 
ycar  period.  SST  analyses  for  ERA-1 5  are  provided  by  the 
United  Kingdom  Meteorological  Office  (UKMO)  for  the 
early  period,  and  by  NOAA  from  November  1981 
onwards.  Sea-ice  cover  has  been  derived  from  satellite  data. 

The  ERA-40  project  applies  a  modem  variational 
data  assimilation  technique  for  the  past  conventional 
and  satellite  observations.  The  model  physics  and  the 
surface  parameterization  have  been  upgraded  and 
improved  since  ERA-15.  A  significant  difference 
between  ERA-40  and  ERA-15  is  in  the  use  of  satellite 
data.  ERA-40  uses  the  Advanced  Tiros  Operational 
Vertical  Sounde  (ATOVS)  radiances  directly,  while  in 
ERA-15  temperature  and  humidity  retrievals  were  used. 
In  addition.  Special  Sensor  Microwave  Imager  (SSM/i) 
and  European  Remote  Sensing  Satellite  (ERS)  data  are 
used  in  ERA-40.  The  SST/lce  data  set  produced  by  the 
Hadley  Centre  and  National  Oceanic  and  Atmospheric 
Administration  (NOAA)ZNational  Environmental  Satel¬ 
lite.  Data,  and  Information  Service  (NESD1S)  has  been 
made  available  to  the  ERA-40  project.  For  use  in 
the  analyses  and  HYCOM  simulations  performed  here, 
25-ycar  monthly  climatologies  of  atmospheric  forcing 
parameters  were  formed  from  the  ERA -40  archives 
spanning  1979-2002  with  the  same  6  hourly  wind 
anomalies  added  as  in  the  ERA-15  ease. 

The  atmospheric  forcing  from  CORE-CNY  is  for  a 
single  year  and  combines  re-analysis  data  from  NCEP 
with  satellite  measurements  to  reduce  errors  existing  in 
the  original  NCEP  fields,  especially  radiation  fields 
(e.g..  Fee  et  at.,  2005).  Heat  fluxes  for  the  CORE-CNY 
forcing  are  produced  by  applying  the  NCAR  bulk 
parameterization  to  NOAA  SSTs,  however  HYCOM 
calculates  its  own  latent  and  sensible  heat  fluxes  based 
on  model  SST  (sec  below).  The  forcing  includes  6 
hourly  representative  variability  in  all  its  fields. 

In  addition  to  wind  and  thermal  forcing,  HYCOM 
forcing  incorporates  monthly  mean  climatologies  of  river 
discharge  and  satellite-based  attenuation  coefficients  for 
Photosynthctically  Active  Radiation  (A>ar  in  m  ').  The 
shortwave  radiation  at  depth  is  calculated  using  a  spatially 
and  temporally  varying  monthly  A>ak  climatology  as 
processed  from  the  daily-averaged  A4W  (attenuation 


coefficient  at  490  nm)  data  set  from  Sca-view  ing  Wide 
Field-of-view  Sensor  (SeaWiFS)  during  1997  2001. 
Thus,  using  ocean  color  data,  the  effects  of  water  turbidity 
are  included  in  the  model  simulations  through  the 
attenuation  depth  (A>ar  in  m)  for  shortwave  radiation. 
The  rate  of  heating/cooling  of  model  layers  in  the  upper 
ocean  is  obtained  from  the  net  heat  flux  absorbed  from  the 
sea  surface  down  to  a  depth,  including  water  turbidity 
effects  (e.g.,  Kara  et  al..  2005a).  Previously,  it  was  shown 
that  water  turbidity  can  be  quite  significant  in  SST 
simulations,  especially  in  equatorial  regions  (e.g.,  Kara 
et  al.,  2004,  2005b). 

The  net  surface  heat  flux  that  has  been  absorbed  (or 
lost)  b>  the  upper  ocean  to  depth  is  parameterized  as  the 
sum  of  the  downward  surface  solar  radiance,  upward 
longwave  radiation,  and  the  downward  latent  and 
sensible  heat  fluxes  (see  Section  2).  Net  solar  radiation 
(the  sum  of  net  shortwave  and  longwave  radiation)  at  the 
sea  surface  is  so  dependent  on  cloudiness  that  it  is  taken 
directly  from  the  given  NWP  product  (ECMWF  or 
CORE-CNY)  for  use  in  the  1 1YCOM.  The  net  longwave 
flux  is  the  sum  of  downward  longwave  (from  the 
atmosphere)  and  upward  black -body  radiation.  The 
NWP  (input)  black-body  radiation  is  corrected  within 
HYCOM  to  allow  for  the  difference  between  NWP  SST, 
7S,  and  HYCOM  SST  (Kara  et  al.,  2(K)5b).  The 
downward  longwave  radiation  is  often  not  archived, 
but  if  archived,  input  net  longwave  is  calculated  using 
an  archived  NWP  SST. 

Latent  and  sensible  heat  fluxes  at  the  air  sea  interface 
are  calculated  using  computationally  efficient  bulk 
formulae  that  include  the  effects  of  dynamic  stability 
(Kara  et  al..  2()()5d).  The  details  of  the  paramctcrizations 
arc  given  in  Section  2.  HYCOM  treats  rivers  as  a  “runoff" 
addition  to  the  surface  precipitation  Held.  The  How  is  first 
applied  to  a  single  ocean  grid  point  and  smoothed  over 
surrounding  ocean  grid  points,  yielding  a  contribution  to 
precipitation  in  m  s  !. 

4.  HYCOM  SST  simulations 

In  this  section,  w'c  investigate  whether  or  not  monthly 
mean  SSI's  obtained  from  atmospherically-forced 
HYCOM,  which  includes  the  effects  of  bulk  parametcr- 
izations  in  its  surface  energy  balance  (see  Section  2),  is 
strongly  controlled  by  the  near-surface  Tm  used  in  the 
atmospheric  forcing.  The  global  HYCOM  simulations 
presented  in  this  paper  were  performed  w  ith  no  assimila¬ 
tion  of  any  oceanic  data  except  initialization  from 
climatology.  1  here  is  only  weak  relaxation  to  sea  surface 
salinity  to  keep  the  evaporation  precipitation  budget  on 
track  in  the  model.  Model  simulations  are  performed 
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using  atmospheric  wind/thermal  forcing  from  each  NWP 
product  (ERA- 1 5,  ERA -40  and  CORE-CN  Y),  separately. 
We  did  not  perform  any  simulations  using  the  alternative 
approach  of  applying  NWP  net  heat  fluxes  and  relaxing  to 
“observed”  SST,  both  because  this  explicitly  includes  the 
“answer”,  at  least  for  SST,  in  the  model  run  and  because 
\vc  would  have  to  choose  one  particular  relaxation  scheme 
which  would  still  leave  open  the  question  of  whether 
some  other  variant  would  be  more  appropriate. 

Near  statistical  equilibrium  was  reached  for  each 
simulation  after  four  model  years.  A  linear  regression 
analysis  was  performed  for  domain  averaged  quantities 
(layer  temperature,  salinity,  potential  and  kinetic  energy, 
etc.)  to  investigate  statistical  equilibrium  in  each  layer, 
and  is  expressed  numerically  as  %  change  per  decade. 
The  model  is  deemed  to  be  in  statistical  equilibrium 
when  the  rate  of  potential  energy  change  is  acceptably 
small  (e.g.,  <  1%  in  5  years)  in  all  layers.  The  statistical 
equilibrium  was  accomplished  using  climatological 
monthly  mean  thermal  atmospheric  forcing,  but  with 
wind  forcing  that  includes  the  6-h  sub-monthly  varia¬ 
bility  because  of  mixed  layer  sensitivity  to  high 
frequency  forcing  (e.g.,  Kara  et  aL  2005a). 

Performing  a  1 -year  global  HYCOM  simulation  using 
any  of  the  atmospheric  forcing  products  required  ~15 
wall-clock  hours  on  64  HP/Compaq  SC45  processors. 
Thus,  the  0.72°  global  HYCOM  provides  inexpensive, 
non-eddy-resolving  ocean  model  simulations. 

Monthly  mean  SST  fields  obtained  from  the  model 
simulations  are  evaluated  through  extensive  model-data 
comparisons  using  various  statistical  metrics  (Section 
4.1 ).  Within  the  quantitative  framework,  statistical  error 
and  skill  analyses  arc  then  presented  for  quality 
assessment  of  the  model  in  simulating  climatological 
monthly  mean  SST  using  bulk  parameterizations  and  the 
three  atmospheric  forcing  products  (Section  4.2). 


are  mainly  designed  for  large-scale  climate  studies. 
They  arc  derived  by  a  linear  interpolation  of  the  w  eekly 
optimal  interpolation  (01),  and  use  in-situ  and  satellite 
SSTs,  making  it  a  reliable  candidate  for  model-data 
comparisons  over  the  global  ocean.  The  existence  of  the 
ice  field  in  the  NOAA  data  set  is  also  an  advantage  for 
the  model  SST  validation  in  the  Arctic  and  Antarctic. 
We  will  also  use  7^  climatologies  from  NWP  products  to 
compare  with  HYCOM  SSTs  in  Section  4.2. 

Monthly  mean  HYCOM  SST  obtained  from 
HYCOM  simulations  using  wind  and  thermal  forcing 
from  ERA-15,  ERA-40  and  CORE-CN Y  is  validated 
using  mean  error  (ME),  root-mean-square  (RMS)  SST 
difference,  correlation  coefficient  (R)  and  non-dimen¬ 
sional  skill  score  (SS).  Let  A',  (C  1 , 2,-,  12)  be  the  set  of 
monthly  mean  NOAA  reference  (observed)  SST  values 
from  January'  to  December,  and  let  Y,  (i~  I,  2,-,  12)  be 
the  set  of  corresponding  HYCOM  estimates  at  a  model 
grid  point.  Also  let  A'  (F)  and  crx  {a })  be  the  mean  and 
standard  deviations  of  the  reference  (estimate)  values, 
respectively.  The  statistical  relationships  (e.g..  Murphy, 
1995)  between  NOAA  and  HYCOM  SST  time  series  at 
a  given  grid  point  are  expressed  as  follows: 

ME -7  A7.  (6) 

1  /2 
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4  1.  Statistical  metrics 

Monthly  mean  HYCOM  SST  climatologies  are 
constructed  from  SSTs  obtained  from  the  last  4  (out  of 
8)  model  years  For  example,  the  climatological  mean 
January  SST  is  formed  using  mean  January  SSTs  from 
model  years  5  8.  The  same  process  is  repeated  for  the 
other  months.  Because  all  simulations  performed  with 
the  0.72°  model  use  climatological  mean  forcing  from 
NWP  products  and  there  are  no  mid-latitude  mesoscalc 
eddies,  the  4-ycar  averaging  period  is  sufficient. 

For  model-data  comparisons  the  NOAA  SSI'  clima¬ 
tology  (Reynolds  et  ah,  2002)  is  taken  as  a  reference 
(truth)  because  its  resolution  (1°*  1°)  is  close  to  that  of 
I IYCOM  (0.72°  x  0.72°  cos(lat)).  The  NOAA  SST  fields 


ME  (i.c.,  annual  bias)  is  the  mean  error  between  the 
mean  HYCOM  and  NOAA  SST  values,  RMS  (root- 
mean-squarc)  SST  difference  is  an  absolute  measure  of 
the  distance  between  the  two  time  scries,  and  the  R 
value  is  a  measure  of  the  degree  of  linear  association 
between  the  time  series.  The  non-dimensional  SS  given 
in  Eq.  (9)  is  the  fraction  of  variance  explained  by 
HYCOM  minus  two  dimensionless  biases  (conditional 
bias,  #c£Mn |,  and  unconditional  bias,  /iulKtind)  which  arc 
not  taken  into  account  in  the  R  formulation  (8)  as 
explained  in  Murphy  (1988),  in  detail.  described 

as  systematic  bias,  is  a  measure  of  the  difference  between 
the  means  of  NOAA  and  HYCOM  SST  time  series.  IK| 
is  a  measure  of  the  relative  amplitude  of  the  variability  in 
the  NOAA  and  1  IYCOM  time  series  or  simply  a  bias  due 
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to  dilTcrcnccs  in  standard  deviations  of  the  SST  time 
series.  Note  SS  is  1.0  (negative)  for  perfect  (poor) 
I  IYCOM  SST  simulations. 

4.2.  IIYCOM  SST  evaluation 

A  model  validation  study  with  respect  to  NOAA  SST 
and  NWP  Ta  from  each  NWP  product  is  presented  using 
the  statistical  metrics  explained  in  Section  4  1.  Our 
puipose  is  to  determine  whether  or  not  HYCOM  SST 


(a)  Bias:  IIYCOM  SST  -  NOAA  SST 


compares  with  the  NOAA  SST  climatology  better  than 
T,  climatologies  from  each  NWP  product.  I  he  original 
monthly  mean  NOAA  SST  climatology  (1°)  was 
interpolated  to  the  global  HYCOM  grid  (0.72°)  for 
comparisons  with  the  model  SSI's.  Similarly,  near- 
surface  ra  from  NWP  products  (ERA-15,  ERA-40  and 
CORE-CNY)  is  also  interpolated  to  the  model  grid. 

The  mean  HYCOM  SST  error,  with  respect  to  the 
NOAA  SST  climatology  (Fig.  4a)  and  Ta  climatologies 
(Fig.  4b)  presents  striking  differences.  In  general, 


<b)  Bias:  IIYCOM  SST  -  NWP  Ta 


big.  4  (a)  Annual  mean  error  (bias)  obtained  from  IIYCOM  simulations  using  climatological  atmospheric  forcing  from  ERA-15,  ERA-40  and 
CORE-CNY,  when  lhe  model  SST  is  compared  lo  (a)  NOAA  SST  climatology  and  (b)  ncar-snrfacc  atr  temperature  climatology  from  each  NWP 
product  Atmospherical  Iv-forccd  IIYCOM  simulations  include  no  assimilation  of  any  SST  or  air  temperature.  In  the  model  simulations,  there  is  no 
relaxation  to  any  SST  or  air  temperature  climatology.  We  evaluate  lime  series  of  HYCOM  SST  versus  NOAA  SST  (similarly  IIYCOM  SST  versus 
NWP  /«)  climatologies  from  January  to  December  at  each  model  grid  point,  producing  lhe  annual  mean  bias  maps  over  the  global  ocean.  In  the  maps, 
the  white  color  represents  bias  values  between  -0.25  CC  and  0,25  °C. 
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most  of  the  mean  errors  shown  in  Fig.  4  do  not  look  like 
the  Ta  7’s  fields,  where  Ts  is  the  SST  used  by  the 
NWP  produets,  presented  in  f  ig.  1  and  discussed  in 
Seetion  2.2. 

To  better  visualize  whether  or  not  HYCOM  mean  SST 
error  patterns  resemble  Tn  SST  patterns  for  the  given  re¬ 
analysis  product,  we  calculate  zonal  averages  of  both 
variables  (F  ig.  5).  Zonaliy-averagcd  Ta  SST  values  can 
also  be  quite  different  for  each  of  these  re -analysis  products, 
it  should  be  emphasized  that  HYCOM  simulations  using 
ERA-15,  ERA-40  and  CORE-CNY  forcing  have  different 
bias  values,  being  wannest  (Table  2)  with  a  globally 
averaged  bias  (HYCOM-NOAA  SST)  value  of().57  °C  for 
the  CORE-CNY  forced  simulation.  In  general,  the  annual 
mean  SST  bias  between  HYCOM  and  the  NOAA 
climatology  is  small  (<0.5  °C)  for  all  three  simulations 
using  the  different  atmospheric  forcing  products  over  most 
of  global  ocean. 

Most  importantly,  there  is  no  clear  relationship 
between  HYCOM  SST  bias  (i.e.,  HYCOM  SST-NOAA 
SST)  and  TiX  SST  fields  (I  ig.  6)  when  using  any  of  the 
atmospheric  re-analysis  products  for  forcing  the  ocean 


model.  This  is  true  despite  the  fact  that  the  HYCOM  SST 
bias  and  7!, -SST  values  from  a  particular  re-analysis 
product  (e.g.,  ERA-15)  ean  be  quite  different  than  those 
from  the  other  product  (e.g.,  ERA -40  and  CORE-CNY). 

If  the  model  SST  is  constrained  by  7j,  from  any 
particular  NWP  product,  then  the  differences  in  ncar- 
surfaee  air  temperatures  between  two  NWP  products 
should  be  similar  to  differences  in  model  SSTs,  obtained 
using  the  corresponding  NWP  forcing  in  simulations. 
This  is  certainly  not  the  case.  As  seen  from  Fig.  7,  the 
scatter  diagram  of  differences  in  T.x  between  ERA-40 
and  CORE-CNY  versus  those  in  SST  from  a  HYCOM 
simulation  forced  with  ERA-40  and  CORE-CNY  are 
quite  different.  While  differences  in  Tn  are  generally 
between  0  °C  and  I  °C,  differences  in  SST  are  generally 
between  1.5°  and  I  °C.  There  is  only  a  weak 
correlation  (/?-~0.36)  between  the  two.  A  similar 
analysis  for  ERA- 1 5  versus  CORE-CNY  is  not  shown 
because  7j,  from  ERA -40  and  ERA- 15  are  almost  same 
over  the  global  ocean. 

As  further  evidence  that  bulk  parametcrizations  are 
not  excessively  tracking  7^,  the  skill  of  monthly  mean 


—  HYCOM  SST  NOAA  SST 

—  NWP  /  HYCOM  SST 


I  .anilide  hell 


Fig  5.  Zonal  averages  of  mean  error  between  HYCOM  anil  the  NOAA  SST  climatology  when  the  model  uses  climatological  atmospheric  forcing 
(wind  and  thermal  forcing)  from  the  three  NW  P  products:  FRA- 1 5,  PRA-40  and  CORF-C'NY.  Also  included  are  zonal  averages  of  Ta  SSI  for  each 
product.  Zonal  averages  arc  calculated  al  each  1°  latitude  hell  over  the  global  ocean 
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Table  2 


Global  av  erages  of  statistical  metrics  calculated  over  the  seasonal  cycle 
between  IIYCOM  and  NOAA  SST  (see  tevl  for  details) 


FRA- 15 

F.RA-40 

CORE-CNY 

NOAA  SST 

Mt  <°C) 

0.23 

0.42 

0.57 

RMS  ( C) 

0  SI 

0.82 

0.98 

^enm  1 

o.ox 

0.07 

0.07 

^tincoml 

0.20 

0.20 

0.42 

R 

0.02 

0.02 

0.92 

SS 

0  48 

040 

0.36 

YH/J  Ta 

MKfC) 

1.54 

1  73 

2.32 

RMS  (CC) 

1.04 

2.08 

2.76 

H 

COlkJ 

0.05 

0.05 

0.07 

^Ciru'oihl 

0.46 

0  54 

0  73 

R 

0.87 

0.86 

0.81 

SS 

0.25 

0.15 

0  14 

Results  arc  given  for  I IYCOM  simulations  using  atmospheric  forcing 
(wind  and  thermal  forcing)  from  HRA-15.  \  RA-40  and  CORE-CNY. 
Basin-averaged  means  are  calculated  over  the  entire  ice-free  global 
IIYCOM  domain.  Similar  comparisons  are  made  between  IIYCOM 
SST  and  NVVP  air  temperature  (a  different  7!,  for  each  product).  A  SS 
value  of  1  indicates  perfect  SST  match  with  respect  to  NOAA  SST  (or 
NWP  Ci)  climatology.  All  R  values  are  statistically  significant  in 
comparison  to  0  70  at  a  95%  confidence  interval. 

model  SSI  against  monthly  NOAA  SST  climatology 
and  against  the  monthly  mean  re -analysis  Ta  (a  different 
Ta  for  caeh  simulation)  is  calculated  (Fig.  8).  We 
evaluate  time  series  of  HYCOM  SST  versus  NOAA 
SST  (similarly  HYCOM  SST  versus  NWP  air  tempera¬ 
ture)  climatologies  from  January  to  December  at  eaeh 
model  grid  point  (see  Section  4.1).  Thus  wc  produce 
spatial  variation  of  the  skill  score  over  the  global  ocean. 
In  the  panels,  blue  (red)  eolor  shows  negative  (positive) 
skill  scores,  and  white  eolor  denotes  skill  score  values 
close  to  0.  In  general,  the  success  of  HYCOM  to  predict 
monthly  SST  on  climatological  time  seales  is  evident 
when  using  atmospheric  forcing  from  any  one  of  the 
NWP  products,  in  that  there  is  positive  skill  over  most  of 
the  glohal  ocean.  The  SST  skill  is  even  close  to  perfect 
(i.e.,  I)  everywhere  except  the  equatorial  Pacific  and 
high  southern  latitudes.  The  low  skill  is  due  in  part  to  the 
atmospheric  forcing  and  in  part  to  deficiencies  in  the 
model,  including  relatively  coarse  model  resolution. 

The  skill  is  much  higher  against  observed  SST  than 
against  the  applied  TA  over  the  most  of  global  ocean 
(fig.  8b).  This  result  is  especially  evident  in  the  Indian 
Ocean  and  Atlantic  Ocean  for  the  simulations  using  the 
different  atmospheric  forcing  products.  When  IIYCOM 
SST  is  validated  against  NOAA  SST  (NWP  Tn\  global 
av  erage  of  SS  values  are  0.48  (0.25),  0.49  (0.15)  and 
0.36  (  0. 1 4 )  for  the  ERA- 1 5,  ERA-40  and  CORE-CN  Y 


forcing  eases,  respectively  ( fable  2).  The  reduction  in 
SS  values  is  generally  >50%  when  HYCOM  SST  is 
evaluated  with  respect  to  NPW  Ta  as  opposed  to  NOAA 
SST. 

The  low  skill  between  IIYCOM  SST  and  NWP  /a  is 
due  mainly  to  the  fact  that  Z?unc ontl  is  significantly 
increased  (not  shown),  i.e.,  there  are  large  biases  in 
means  of  HYCOM  SST  and  NWP  based  on  the 
definition  of  unconditional  bias  (Section  4.1).  Note  that 
there  is  not  trtueh  change  in  global  averages  of  R  and 
Bcond.  Further  comparisons  of  HYCOM  SST  versus 
both  NOAA  SST  and  NWP  Ta  elearly  demonstrate  that 
the  shape  of  the  seasonal  cycle  between  two  variables 
follow  eaeh  other  quite  well  as  evident  from  R  values 
>0.8  (Fig.  9).  The  exception  is  the  I  IYCOM  simulation 
with  SST  compared  to  T.A  from  CORE-CNY.  As 
expected,  all  NWP  models  have  different  boundary 
layer  parameterizations.  In  addition,  data  assimilation 
methods  and  data  type  (e.g.,  satellite  and  in-situ)  used  in 
the  assimilations  may  differ  from  one  NWP  center  to 
another  one.  Therefore,  differences  in  their  outputs,  sneh 
as  the  near-surface  air  temperatures  are  expected. 

Finally,  overall  IIYCOM  performance  in  simulating 
cliniatologieal  SST  is  diseussed  when  the  model  uses 
atmospheric  forcing  from  ERA-15,  ERA-40  and  CORE- 
CNY.  Table  2  provides  a  summary  of  statistical  metries 
between  HYCOM  SST  and  NOAA  SST  in  the  form  of 
global  averages.  While  the  SST  simulated  by  HYCOM 
using  atmospheric  forcing  from  CORE-CNY  did  not  yield 
results  as  good  as  those  using  ERA- 1 5  and  ERA-40, 
HYCOM  generally  gave  similar  sueecss  when  using  any 
of  the  atmospheric  forcing  produets. 

Although  there  have  been  many  improvements  in  the 
ERA -40  forcing  in  comparison  to  the  I  RA- 15  forcing 
(Section  3.2),  the  model  response  in  simulating 
ehmatologieal  SST  remains  similar  with  global  mean 
RMS  SST  differences  of  0.81  °C  and  0.82  °C,  and  SS 
values  of  0.48  and  0.49  for  the  ERA- 15  and  ERA-40 
cases.  There  is  a  warm  model  SST  bias  in  all  of  the 
simulations.  This  could  he  due  to  known  shortcomings 
in  NWP  solar  radiation  fields,  e.g.,  due  to  inaccurate 
cloudiness,  or  could  result  from  the  ocean  model 
climatology.  In  any  case,  such  a  bias  can  be  removed 
by  applying  a  spatially  varying  but  constant  in  time 
correction  to  the  total  heal  fluxes  (not  diseussed  here). 

5.  Conclusions 

The  most  significant  result  of  this  paper  is  that 
virtually  all  applications  of  the  bulk  formulae,  including 
a  fixed  air  temperature,  does  not  make  the  sensible  and 
latent  heat  behave  as  though  the  model  SST  tracks  the 
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l;ig.  6.  Scatter  plots  for  I IYCOV1  SST  bias  (HYC'OM  SST-NOAA  SST)  versus  f  SSI*  based  on  zonal ly-averagcd  values  shown  in  f  ig.  5.  The  least 
square  line  for  each  product  is  also  shown  along  with  the  linear  R  value  between  the  two  variables  for  each  product,  respectively.  TIk  least  squares 
lines  are  based  on  zonal ly-aver«igcd  \  allies  for  each  NWP  product.  Zonally-averaged  values  are  used  for  plotting  purposes  to  reduce  number  of  data 
over  the  global  ocean.  The  regions  where  ice  is  located  are  not  used  in  the  analyses. 
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I  ig.  7.  Scatter  plots  for  SS  I  differences  from  two  HYC’OM  simulations  forced  with  HRA-40  and  C'ORl  -('NY  versus  differences  in  berween  l  RA- 
40  and  CORH-C'Nf  Results  are  hased  on  annual  mean  values  over  the  global  ocean  except  ice-covcred  regions.  The  least  squares  line  having  a  linear 
R  value  of  0.36  is  also  show  n 
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1  ig.  X.  Dimensionless  skill  score  calculated  over  the  seasonal  cycle  for  I IYCOM  SS  I  with  respect  to  (a)  NOAA  SS  I  climatology  and  (b)  near-surface 
air  temperature  climatology  from  each  NWP  product  Results  arc  shown  when  the  model  is  forced  with  NWP  products  (KRA- 1 5.  HRA-40  and  CORK  - 
t  NY).  separately  Note  that  air  temperature  climatology  used  for  comparisons  tn  (b)  is  different  for  each  panel.  As  tn  fig  4.  atmospherically-forced 
1  IYCOM  simulations  include  no  assimilation  of  any  SST  or  air  temperature,  also  no  relaxation  to  any  SST  or  air  temperature  climatology. 


near-surfaee  air  temperature  too  strongly.  This  is 
demonstrated  using  an  atmospherieally-foreed  OGCM 
(1  IYCOM)  with  no  assimilation  of  any  SST  and  no 
relaxation  to  any  SST  climatology. 

Accurate  SST  simulations  in  an  OGCM  depend  on 
how  the  local  changes  in  SST  made  by  the  bulk  formula 
arc  modified  by  vertical  mixing  (mixed  layer  physics) 
and  adveetion.  Atmosphcrically-Iorccd  1  IYCOM  simu¬ 
lations  clearly  reveal  that  while  the  RMS  SS  I  difference 
is  quite  small  between  SST  and  near-surfaee  air 
temperature  over  most  of  the  global  ocean,  the  use  of 


the  latter  in  a  hulk  parameterization  does  not  introduce 
an  inappropriate  restoring  toward  observed  SST  in  the 
model  simulations. 

Through  extensive  model-data  comparisons  using 
various  statistical  metrics,  we  have  demonstrated  that 
bulk  parameterization^  track  observed  SST  rather  than 
near-surfaee  air  temperature  in  atmospherieally-foreed 
OGCM  simulations.  When  monthly  mean  SST  simul¬ 
ated  by  HYCOM  is  compared  against  the  monthly 
mean  observed  SST  climatology  from  NOAA  and 
climatologies  from  ERA-15,  ERA-40  and  CORE-CNY, 
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Tig.  9.  The  same  as  Tig.  S  hut  for  the  correlation  coefficient.  In  the  maps,  the  white  color  represents  correlation  values  of  >0.95 


the  model  skill  is  higher  against  observed  SST  than 
against  the  applied  7j  used  in  the  simulations.  Overall,  the 
global  mean  of  RMS  SST  dilTcrence  For  HYCOM  is 
0.81  °C  0.82  °C  and  0,98  °C  with  respect  to  the  NOAA 
climatology  over  the  seasonal  cycle  For  the  simulations 
using  atmospheric  w  ind  and  thermal  forcing  from  ERA- 
1 5,  ERA-40  and  CORE-CNY,  respectively.  These  global 
mean  RMS  differences  increase  by  ~  150%  (1.94  °C, 
2.08  °C  and  2.76  °C,  respectively)  when  HYCOM  SS  I  is 
evaluated  with  against  near-surfacc  air  temperature  from 
the  three  NWP  produets.  Model  SST  validation  against 
NOAA  SST  and  NWP  7*  yields  similar  correlation 
coefficients  (generally  >0.90)  over  the  most  of  global 
ocean,  w  hich  may  imply  an  indirect  relaxation  of  model 
SST  to  NWP  Ta.  However,  significantly  large  differences 


in  dimensionless  skill  score  and  RMS  differences 
between  the  pairs  of  HYCOM  SST  versus  NOAA  SST 
and  model  SST  versus  NWP  Ta  indicate  that  such  a  direct 
relaxation  is  not  the  case.  Further,  the  correlation 
between  Ta(ERA-40)-ra(CORE-CNY)  and  HYCOM 
SST(ERA-40)-H YCOM  SST(CORE-CNY)  is  R  = 
0.36  when  SST  relaxation  to  NWP  near-surface  air 
temperature  would  imply  a  substantial  positive  value. 
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\ppcndix  A.  New  features  of  HYCOM 

HYCOM  development  has  continued  since  the  first 
release  discussed  in  Bleck  (2002),  and  new  features 
have  been  added  to  the  model.  A  short  summary  of  the 
new  improvements  in  the  model  code  (version  2. 1.18)  is 
provided  here. 

A  hybrid  vertical  coordinate  grid  generator  (called 
hybgcn)  is  used  in  choosing  the  optimal  vertical 
coordinate  at  each  location  every  time  step.  It  then 
remaps  the  vertical  coordinate  accordingly.  I  here  are 
significant  improvements  to  hybgcn  in  the  latest  version 
of  the  model.  The  remapper  in  the  original  model 
assumed  that  each  field  was  constant  in  the  vertical 
within  each  layer.  When  remapping  layers  that  arc  far 
from  the  target  isopyenal,  this  approach  can  lead  to 
excessive  diffusion.  The  current  modified  remapper  in 
HYCOM  allows  the  profile  to  vary  linearly  across  a 
layer  if  the  layer  is  not  close  to  being  isopyenal,  which 
significantly  reduces  diffusion.  In  such  cases,  the 
original  rcmappei  used  donor-cell  upwind  advection, 
and  the  latest  remapper  uses  the  piecewise  linear  method 
(van  Leer  1977;  Lin  et  al.,  1994)  with  a  monotomzed 
centra  I -difference  limiter.  The  hybrid  grid  generator  and 
horizontal  advection  terms  can  each  conserve  potential 
temperature  and  salinity,  potential  density  and  salinity  or 
potential  density  and  potential  temperature.  Because  the 
hybrid  grid  generator  must  act  in  density  space,  potential 
density  and  salinity  arc  topically  conserved  in  both  the 
grid  generator  and  advection. 

Since  the  first  release  of  HYCOM  (Bleck,  2002), 
there  is  improved  support  lorz  and  o  levels  in  shallow 
water,  and  more  flexible  selection  of  Laplacian  and 
Inharmonic  diffusion  There  arc  several  scalar  advection 
techniques  available  in  HYCOM,  such  as  the  donor  cell, 
the  Flux-Corrected  Transport  (FCT)  scheme  (2nd  and 
4th  order)  and  Multidimensional  Positive  Definite 
Advection  Transport  Algorithm  (MPDATA)  (e.g., 
Smolarkiewiez  and  Margolin,  1998;  Balsara  and  Shu. 
2000).  While  we  do  not  discuss  it  here,  HYCOM 
simulations  demonstrate  that  FCT  is  more  computation¬ 


ally  efficient  and  less  diffusive  than  MPDATA, 
consistent  with  results  by  Smolarkiewiez  (1984).  Multi¬ 
ple  tracers  and  off-line  one-way  nesting  are  possible  in 
the  model.  HYCOM  includes  a  hybrid  to  fixed  vertical 
grid  remapper,  allowing  fixed-coordinate  nests  inside 
hybrid  coordinate  outer  domains.  Vertical  remapping 
uses  the  Piecewise  Linear  Method  (PLM)  for  fixed- 
coordinate  layers.  Isopyenal  layer  target  densities  can 
vary  spatially,  and  stability  of  these  layers  is  controlled 
by  locally-referenced  potential  density. 

HYCOM  can  be  run  using  any  one  of  five  mixed 
layer  models:  (I)  The  K-Profile  Parameterization  (KPP) 
level  I  turbulence  closure  (Large  el  al.,  1997),  (2) 
Kraus  Turner  (KT)  mixed  layer  model  (Kraus  and 
Turner,  1967),  (3)  Price,  Weller  and  Pinkel  (PWP) 
mixed  layer  model  (Price  et  al.,  1986),  (4)  the  Mel  lor 
Yaniada  (MY)  level  2.5  turbulence  closure  (Mcllor  and 
Yamada,  1982),  and  (5)  the  NASA  Goddard  Institute  for 
Space  Studies  (GISS)  level  2  turbulence  closure  (Canute 
et  al.,  2002).  I  he  model  can  also  be  run  with  no  mixed 
layer  model. 

Earlier  HYCOM  simulations  used  the  bulk  parameter- 
izations  presented  in  Kara  et  al.  (2002),  which  has  a  few 
shortcomings  in  the  determination  of  exchange  coeffi¬ 
cients  of  sensible  and  latent  heat  fluxes  (Cs  and  C|)  at  high 
and  low  wind  speeds.  This  was  due  to  limitations  in  the 
observation-based  input  data  sets  used  for  deriv  ing  the 
exchange  coefficients.  Also,  the  effect  of  humidity  was 
not  included  in  the  parameterizations  of  air  sea  stability. 
With  the  availability  of  the  improved  COARI  algorithm 
version  3.0  (Fairall  et  al.,  2003),  Kara  et  al.  (2005d) 
derived  new  polynomials  for  C|  and  Cs  for  use  in 
HYCOM,  including  the  full  effect  of  stability  for  a  w  ide 
range  of  conditions  occurring  over  the  global  ocean,  flic 
new  data  set  from  C'OARE  3.0  reduces  the  under¬ 
estimation  or  overestimation  in  the  polynomials  10  to 
20%)  used  for  deriving  the  exchange  coefficients 
presented  in  Kara  et  al  (2000,  2002).  In  addition,  they 
have  the  advantage  of  providing  reliable  transfer 
coefficients  at  low  and  high  wind  speeds.  They  represent 
only  an  approximation  to  the  COARL  algorithm,  but  with 
the  advantage  of  robustness  and  computational  efficiency 
that  make  it  suitable  for  use  in  various  air  sea  interaction 
applications  and  in  any  OGCM. 
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